// A Javascript implementation of the Nelder-Mead Simplex method, mostly copied
// from
//     http://shell.lab49.com/~vivake/python/Simplex.py

function Simplex(optFun, initial, stepSizes){
	this.optFun = optFun;
	this.guess = initial;
	this.stepSizes = stepSizes;
	this.kR = -1;  // reflection constant
	this.kE = 2;   // expansion constant
	this.kC = 0.5; // contraction constant
	this.numvars = initial.length;
	this.simplex = new Array();
	
	this.lowest = -1;
	this.highest = -1;
	this.secondhighest = -1;
	
	this.errors = new Array();
	this.currenterror = 0;
	
	this.done = false;
	// initialize vertices
	for (var vertex = 0; vertex < this.numvars + 3; vertex++){ // two extras to store centroid and reflected point.
		this.simplex[vertex] = new Array();
		for (var x = 0; x < this.numvars; x++){
			this.simplex[vertex][x] = this.guess[x];
		}
	}
	for (var vertex = 0; vertex < this.numvars + 1; vertex++){
		for (var x = 0; x < this.numvars; x++){
			if (x == (vertex - 1)){
				this.simplex[vertex][x] = this.guess[x] + this.stepSizes[x];
			}
		}
		this.errors[vertex] = 0;
	}
	this.calculate_errors_at_vertices();
	this.done = false;
	
}

Simplex.prototype.calculate_errors_at_vertices = function(){
	for (var vertex = 0; vertex < this.numvars + 1; vertex++){
		if (vertex == this.lowest){
			continue;
		}
		for (var x = 0; x < this.numvars; x++){
			this.guess[x] = this.simplex[vertex][x];
		}
		this.currenterror = this.optFun(this.guess);
		this.errors[vertex] = this.currenterror;
	}
}

Simplex.prototype.minimize = function(epsilon, maxiters, monitor){
	this.done = false;
	// walks the simplex down to a local minima.
	for (var iter = 0; iter < maxiters; iter++){
		if (monitor){
			document.getElementById('simplex').innerHTML = document.getElementById('simplex').innerHTML + iter + ':'
		}
		
		
		this.done = this.minimizeIter(epsilon, monitor);
		if (this.done){
			if (monitor){
				document.getElementById('simplex').innerHTML = document.getElementById('simplex').innerHTML +
			          iter + ' iterations, optimum of '  + this.currenterror + ' at ' + this.guess;
			}
			
			break;
		}
	}
	for (var x = 0; x < this.numvars; x++){
		this.guess[x] = this.simplex[this.lowest][x];
	}
	this.currenterror = this.errors[this.lowest];
	return this.guess;
}

Simplex.prototype.minimizeIter = function(epsilon, monitor){
	// does one iteration of the Nelder-Mead simplex algorithm.
	this.highest = 0;
	this.lowest = 0;
	for (var vertex = 0; vertex < this.numvars + 1; vertex++){
		if (this.errors[vertex] > this.errors[this.highest]){
			this.highest = vertex;
		}
		if (this.errors[vertex] < this.errors[this.lowest]){
			this.lowest = vertex;
		}
	}
	this.secondhighest = 0;
	for (var vertex = 0; vertex < this.numvars + 1; vertex++){
		if (vertex == this.highest){
			continue;
		}
		if (this.errors[vertex] > this.errors[this.secondhighest]){
			this.secondhighest = vertex;
		}
	}
	
	// test for convergence.
	var S = 0.0;
    var S1 = 0.0;
	
	for (var vertex = 0; vertex < this.numvars + 1; vertex++){
		S = S + this.errors[vertex];
	}
	var F2 = S / (this.numvars + 1);
	for (var vertex = 0; vertex < this.numvars + 1; vertex++){
		var sq = (this.errors[vertex] - F2)
		S1 = S1 + (sq * sq);
	}
	var T = Math.sqrt(S1 / this.numvars);
	
	if (monitor){
		document.getElementById('simplex').innerHTML = document.getElementById('simplex').innerHTML + 
			'Best = ' + this.errors[this.lowest] + ', Worst = ' + this.errors[this.highest] + '<br/>';
	}
	
	if (T <= epsilon){
		// convergence!
		if (monitor){
			document.getElementById('simplex').innerHTML = document.getElementById('simplex').innerHTML + 
			'converged. <br/>';
		}
		return true;
	}
	// didn't converge.
	for (var x = 0; x < this.numvars; x++){
		var S = 0.0;
		for (var vertex = 0; vertex < this.numvars + 1; vertex++){
			if (vertex == this.highest){
				continue;
			}
			S = S + this.simplex[vertex][x];
		}
		this.simplex[this.numvars+1][x] = S / this.numvars;
	}
	
	this.reflect_simplex();
	this.currenterror = this.optFun(this.guess);
	
	if (this.currenterror < this.errors[this.lowest]){
		var tmp = this.currenterror;
		this.expand_simplex();
		this.currenterror = this.optFun(this.guess);
		if (this.currenterror < tmp){
			this.accept_expanded_point()
		} else {
			this.currenterror = tmp;
			this.accept_reflected_point();
		}
	} else if (this.currenterror <= this.errors[this.secondhighest]){
		this.accept_reflected_point();
	} else if (this.currenterror <= this.errors[this.highest]){
		this.accept_reflected_point();
		
		this.contract_simplex();
		this.currenterror = this.optFun(this.guess);
		if (this.currenterror < this.errors[this.highest]){
			this.accept_contracted_point();
		} else {
			this.multiple_contract_simplex();
		}
	} else if (this.currenterror >= this.errors[this.highest]){
		this.contract_simplex()
		this.currenterror = this.optFun(this.guess)
		if (this.currenterror < this.errors[this.highest]){
			this.accept_contracted_point()
		} else {
			this.multiple_contract_simplex()
		}
	}
	return false;
}

Simplex.prototype.contract_simplex = function(){
	for (var x = 0; x < this.numvars; x++){
		this.guess[x] = this.kC * this.simplex[this.highest][x] + (1 - this.kC) * this.simplex[this.numvars + 1][x];
	}
}

Simplex.prototype.expand_simplex = function(){
	for (var x = 0; x < this.numvars; x++){
		this.guess[x] = this.kE * this.guess[x]                 + (1 - this.kE) * this.simplex[this.numvars + 1][x]
	}
}

Simplex.prototype.reflect_simplex = function(){
	for (var x = 0; x < this.numvars; x++){
		this.guess[x] = this.kR * this.simplex[this.highest][x] + (1 - this.kR) * this.simplex[this.numvars + 1][x]
		this.simplex[this.numvars + 2][x] = this.guess[x]  // remember the reflected point.
	}
}

Simplex.prototype.multiple_contract_simplex = function(){
	for (var vertex = 0; vertex < this.numvars + 1; vertex++){
		if (vertex == this.lowest){
			continue;
		}
		for (var x = 0; x < this.numvars; x++){
			this.simplex[vertex][x] = 0.5 * (this.simplex[vertex][x] + this.simplex[this.lowest][x]);
		}
	}
	this.calculate_errors_at_vertices();
}

Simplex.prototype.accept_contracted_point = function(){
	this.errors[this.highest] = this.currenterror;
	for (var x = 0; x < this.numvars; x++){
		this.simplex[this.highest][x] = this.guess[x];
	}
}

Simplex.prototype.accept_expanded_point = function(){
	this.errors[this.highest] = this.currenterror;
	for (var x = 0; x < this.numvars; x++){
		this.simplex[this.highest][x] = this.guess[x];
	}
}

Simplex.prototype.accept_reflected_point = function(){
	this.errors[this.highest] = this.currenterror;
	for (var x = 0; x < this.numvars; x++){
		this.simplex[this.highest][x] = this.simplex[this.numvars + 2][x];
	}
}